#include <iostream>
#include <algorithm>
#include <vector>
#include <array>
#include <numeric>
#include <functional>
#include "euler/modular.hpp"
#include "euler.h"

BEGIN_PROBLEM(271, solve_problem_271)
	PROBLEM_TITLE("Modular Cubes, Part 1")
	PROBLEM_ANSWER("4617456485273129588")
	PROBLEM_DIFFICULTY(2)
	PROBLEM_FUN_LEVEL(1)
	PROBLEM_TIME_COMPLEXITY("M*3^M")
	PROBLEM_SPACE_COMPLEXITY("M")
	PROBLEM_KEYWORDS("congruence")
END_PROBLEM()

static void solve_problem_271()
{
#if 0
	const int M = 2;
	const std::array<int,M> factors = {{ 7, 13 }};
#else
	const int M = 14;
	const std::array<int,M> factors = {{ 2,3,5,7,11,13,17,19,23,29,31,37,41,43 }};
#endif

	bool verbose = false;
	long long N = std::accumulate(factors.cbegin(), factors.cend(), 
		1LL, std::multiplies<long long>());

	// Solve (x-1)(x^2+x+1) = 0 (mod p) for each prime factor p.
	std::vector<std::vector<int>> roots(M);
	for (int i = 0; i < M; i++)
	{
		int p = factors[i];
		roots[i].push_back(1);
		for (int x = 2; x <= (p - 1) / 2; x++)
		{
			if ((x*x+x+1) % p == 0)
			{
				roots[i].push_back(x);
				roots[i].push_back(p-1-x);
				break;
			}
		}
		if (verbose)
		{
			std::cout << "Cubic roots of 1 mod " << p << " are:";
			std::for_each(roots[i].cbegin(), roots[i].cend(), [](int x) { 
				std::cout << " " << x;
			});
			std::cout << std::endl;
		}
	}

	// Pre-compute the following for Chinese remainder theorem.
	// b = modinv(N/p_i, p_i)
	std::vector<long long> b(M);
	for (int i = 0; i < M; i++)
		b[i] = euler::modinv(N / factors[i], (long long)factors[i]);

	// For each combination of a_i, i=1..m, solve the linear congruence 
	// system: x = a_i (mod p_i).
	std::vector<int> choice(M);
	long long sum = 0;
	while (true)
	{
		// Compute the root.
		long long x = 0;
		for (int i = 0; i < M; i++)
		{
			long long t = roots[i][choice[i]] * b[i] * (N / factors[i]);
			x += t;
			x %= N;
		}
		sum += x;
		if (verbose)
			std::cout << "Root: " << x << std::endl;

		// Get next permutation
		int i;
		for (i = 0; i < M; i++)
		{
			if (++choice[i] < (int)roots[i].size())
				break;
			choice[i] = 0;
		}
		if (i >= M)
			break;
	}
	std::cout << (sum - 1) << std::endl;
}
